Efficient posterior profiles

Marco Raveri (marco.raveri@unige.it), Cyrille Doux (doux@lpsc.in2p3.fr), Shivam Pandey (shivampcosmo@gmail.com)

In this notebook we show how to obtain posterior profiles from synthetic probability models, as in Raveri, Doux and Pandey (2024), arXiv:2409.09101.

If you want more details on how to build normalizing flow based synthetic models for posterior distributions check out the corresponding example notebook.

Table of contents

1. Notebook setup
2. Flow training
3. Posterior profiles
4. Profile accuracy
5. Real world application
   1. Best constrained parameters profile
   2. Full profile triangle plot

Notebook setup:

We start by importing everything and setting up a controlled example:

# Show plots inline, and load main getdist plot module and samples class
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext autoreload
%autoreload 2

# import libraries:
import sys, os
os.environ['TF_USE_LEGACY_KERAS'] = '1'  # needed for tensorflow KERAS compatibility
os.environ['DISPLAY'] = 'inline'  # hack to get getdist working
sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'../..')))
from getdist import plots, MCSamples
from getdist.gaussian_mixtures import GaussianND
import getdist
getdist.chains.print_load_details = False
import scipy
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# tensorflow imports:
import tensorflow as tf
import tensorflow_probability as tfp

# import the tensiometer tools that we need:
import tensiometer
from tensiometer.utilities import stats_utilities as utilities
from tensiometer.synthetic_probability import synthetic_probability as synprob

# getdist settings to ensure consistency of plots:
getdist_settings = {'ignore_rows': 0.0, 
                    'smooth_scale_2D': 0.3,
                    'smooth_scale_1D': 0.3,
                    }    

2024-12-11 13:20:02.390337: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.

We build here a random Gaussian mixture model that we are going to use for tests.

The esample is seeded so that it is reproducible. If you want a different example change the value of the seed.

You can also change dimensionality and number of modes.

# define the parameters of the problem:
dim = 6
num_gaussians = 3
num_samples = 10000

# we seed the random number generator to get reproducible results:
seed = 100
np.random.seed(seed)
# we define the range for the means and covariances:
mean_range = (-0.5, 0.5)
cov_scale = 0.4**2
# generate means and covs:
means = np.random.uniform(mean_range[0], mean_range[1], num_gaussians*dim).reshape(num_gaussians, dim)
weights = np.random.rand(num_gaussians)
weights = weights / np.sum(weights)
covs = [cov_scale*utilities.vector_to_PDM(np.random.rand(int(dim*(dim+1)/2))) for _ in range(num_gaussians)]

# cast to required precision:
means = means.astype(np.float32)
weights = weights.astype(np.float32)
covs = [cov.astype(np.float32) for cov in covs]

# initialize distribution:
distribution = tfp.distributions.Mixture(
    cat=tfp.distributions.Categorical(probs=weights),
    components=[
        tfp.distributions.MultivariateNormalTriL(loc=_m, scale_tril=tf.linalg.cholesky(_c))
        for _m, _c in zip(means, covs)
    ], 
    name='Mixture'
    )

# sample the distribution:
samples = distribution.sample(num_samples).numpy()
# calculate log posteriors:
logP = distribution.log_prob(samples).numpy()

# create MCSamples from the samples:
chain = MCSamples(samples=samples, 
                  settings=getdist_settings,
                  loglikes=-logP,
                  name_tag='Mixture',
                  )

# we want to find the maximum posterior point:
_temp_maxima, _temp_max_value = [], []
for _m in means:
    # maximize the likelihood starting from all the means:
    _max = scipy.optimize.minimize(lambda x: -distribution.log_prob(x).numpy(), _m, method='Nelder-Mead')
    # this usually converges to the nearest mode:
    _temp_maxima.append(_max.x)
    _temp_max_value.append(-_max.fun)
maximum = _temp_maxima[np.argmax(_temp_max_value)]
maximum_value = _temp_max_value[np.argmax(_temp_max_value)]

# we make a sanity check plot:
g = plots.get_subplot_plotter()
g.triangle_plot(chain, filled=True, markers=maximum)

As we can see this example beautifully showcases projection effects, especially in $p_5$. As you might have seen in the synthetic probability notebook the low peak in the 1D marginal of $p_5$ is the actual full-D peak.

Flow training:

We now train a synthetic probability model. Note that we need a flow with good local accuracy since we are going to maximize its value.

If you are interested in how to build and control these types of flows, check out the synthetic probability tutorial.

kwargs = {
          'feedback': 1,
          'verbose': -1,
          'plot_every': 1000,
          'pop_size': 1,
          'num_flows': 5,
          'epochs': 300,
        }

flow = synprob.average_flow_from_chain(chain,  # parameter difference chain
                                       **kwargs)

Training flow 0

0epoch [00:00, ?epoch/s]

Training flow 1

0epoch [00:00, ?epoch/s]

Training flow 2

0epoch [00:00, ?epoch/s]

Training flow 3

0epoch [00:00, ?epoch/s]

Training flow 4

0epoch [00:00, ?epoch/s]

# we need to check metrics to make sure everything looks good:
flow.print_training_summary()

Number of flows: 5
Flow weights   : [0.2 0.2 0.2 0.2 0.2]
loss         : [7.29 7.31 7.34 7.34 7.38]
val_loss     : [7.31 7.34 7.32 7.33 7.35]
lr           : [3.16e-05 1.00e-03 1.00e-03 3.16e-05 3.16e-04]
chi2Z_ks     : [0.05 0.03 0.03 0.04 0.05]
chi2Z_ks_p   : [0.03 0.31 0.19 0.09 0.03]
loss_rate    : [-1.95e-04 -3.59e-03 -2.08e-03  9.54e-06 -6.07e-04]
val_loss_rate: [ 0.    0.01 -0.01 -0.   -0.  ]

# then we need to check the estimate of the local accuracy:
ev, eer = flow.evidence()
print(f'log(Z) = {ev} +- {eer}')

log(Z) = 0.09110794216394424 +- 0.5065720677375793

It is usually good to check that the evidence error (if available) is under control, especially when we want to obtain posterior profiles.

The evidence error is the variance of log likelihood values on samples from the flow and gives us a handle on the local accuracy of the flow.

Cathastrofic initialization of weights happens and if this value is too high then it might be worth re-running flow training.

If you are training on marginals (without likelihood values available) then it might be a good idea to add population selection as a layer of protection against bad initialization.

If you find that the results are not stable (especially in the bulk of the posterior) check out the notebook on synthetic probability modelling and the section discussing high accuracy flows.

Posterior profile:

We now want to calculate posterior profiles for our distribution. These are obtained maximizing over all parameters but the ones that are been considered.

Having efficient flow models from which we can sample and calculate probability values means that we can afford lots of maximizations, so that we can calculate up to 2D profiles.

Note that a flow can be trained on marginal distributions, allowing us to combine profiling and marginalization.

Posterior profiles can be easily obtained using the appropriate tensiometer class operating on the flow:

from tensiometer.synthetic_probability import flow_profiler

# define the options for the profiler:
profiler_options = {
        'num_gd_interactions_1D': 100,  # number of gradient descent interactions for the 1D profile
        'num_gd_interactions_2D': 100,  # number of gradient descent interactions for the 2D profile
        'scipy_options': {  # options for the scipy polishing minimizer
                    'ftol': 1.e-06,
                    'gtol': 0.0,
                    'maxls': 40,
                },
        'scipy_use_jac': True,  # use the jacobian in the minimizer
        'num_points_1D': 64, # number of points for the 1D profile
        'num_points_2D': 32, # number of points per dimension for the 2D profile
        'smooth_scale_1D': 0.2, # smoothing scale for the 1D profile
        'smooth_scale_2D': 0.2, # smoothing scale for the 2D profile
        }

# initialize the profiler:
flow_profile = flow_profiler.posterior_profile_plotter(flow, 
                                                       initialize_cache=False,  # usually we want to avoid initializing the cache for all the parameters
                                                       feedback=2,  # we want high feedback to see the progress
                                                )

Flow profiler
  * use_scipy  = True
  * pre_polish = True
  * polish     = True
  * smoothing  = True
  * box_prior  = False
  * parameters = ['param1', 'param2', 'param3', 'param4', 'param5', 'param6']
  * feedback   = 2

We now have initialized the profiler. If cache initialization is enabled the code will calculate 1D and 2D profiles for all parameter combinations. This is usually a lot of work so we keep it separate in this tutorial and proceed with the profile calculation for just a few parameters:

# now we initialize the cache for the parameters we want to profile:
profile_params = ['param3', 'param5']
flow_profile.update_cache(params=profile_params, 
                          **profiler_options)

  * initializing profiler data.
    * generating samples for profiler
    - number of random search samples = 20000
    - sampling the distribution
    - time taken to sample the distribution 7.856 (s)
  * finding MAP
    * finding global best fit
    - doing initial randomized search
    - time taken for random initial selection 0.001371 (s)
    - doing minimization (scipy)
    - time taken for minimization 2.085 (s)
  * initializing 1D profiles
    1D profiles:   0%|          | 0/2 [00:00<?, ?it/s]
    * calculating the 1D profile for: param3 with index 2
    - doing initial randomized search
    - time taken for random algorithm 0.9823 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 3 / 61
    - doing gradient descent pre-polishing
    - time taken for gradient descents 35.15 (s)
      at the end of descent after 100 steps 7 samples were still beeing optimized.
      gradient descents did not improve 30 points
    - doing minimization polishing (scipy)
    - time taken for polishing 7.392 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.003166 (s)
    - smoothing results
    - time taken for smoothing 0.001274 (s)
    1D profiles:  50%|█████     | 1/2 [00:43<00:43, 43.53s/it]
    * calculating the 1D profile for: param5 with index 4
    - doing initial randomized search
    - time taken for random algorithm 0.002331 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 2 / 62
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.203 (s)
      at the end of descent after 100 steps 16 samples were still beeing optimized.
      gradient descents did not improve 37 points
    - doing minimization polishing (scipy)
    - time taken for polishing 9.693 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.000474 (s)
    - smoothing results
    - time taken for smoothing 0.001191 (s)
    1D profiles: 100%|██████████| 2/2 [00:54<00:00, 27.22s/it]
  * initializing 2D profiles
    2D profiles:   0%|          | 0/1 [00:00<?, ?it/s]
    * calculating the 2D profile for: param3, param5 with indexes  2 4
    - doing initial randomized search
    - time taken for random algorithm 0.5556 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 385 / 639
    - doing gradient descent pre-polishing
    - time taken for gradient descents 2.386 (s)
      at the end of descent after 100 steps 118 samples were still beeing optimized.
      gradient descents did not improve 203 points
    - doing minimization polishing (scipy)
    - time taken for polishing 72.03 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0181 (s)
    - smoothing results
    - time taken for smoothing 0.002401 (s)
    2D profiles: 100%|██████████| 1/1 [01:14<00:00, 75.00s/it]

As you can notice the profile calculation is complicated and intensive. For this reason caching is implemented and thoroughly used.

After calculating profiles the result can be saved to file for effective caching:

# this command will save the profile to a pickle file:
#flow_profile.savePickle('flow_profile.pkl')
# note that the flow cannot be pickled easily and has its own save and load functions. This means you have to save it separately.
#flow_profile = flow_profiler.posterior_profile_plotter.loadPickle('flow_profile.pkl', flow)

The profiler class hijacks getdist MCSamples so that it can be directly used for getdist plotting as follows:

g = plots.get_subplot_plotter()
g.triangle_plot([flow_profile, flow.MCSamples(10000)], 
                params=profile_params,
                markers=[flow_profile.flow_MAP[flow_profile.index[_p]] for _p in profile_params],
                filled=False, 
                shaded=True, 
                diag1d_kwargs={'normalized':True},
                legend_labels=['Profile','Marginal'])

As we can see the projection effect in $p_5$ is rightfully recovered and the mode that is smaller in the marginal is much higher in the profile.

Note that there might be some little discrepancy between the peak of the profiles and the full-D peak, while there should be no difference. This is due to the finite resolution of the 1D and 2D profiles, which are typically only computed on a small-ish grid.

The profiler class is fully interfaced with getdist plotting facilities. Everything that is not previously cached is recomputed on the flight, as we can see in the following example:

plot_params = ['param1'] + profile_params 
g = plots.get_subplot_plotter()
g.plots_1d([flow_profile, flow.MCSamples(10000)], 
           params=plot_params, 
           legend_labels=['Profile','Marginal'], 
           nx=3, normalized=True,
           markers=[flow_profile.flow_MAP[flow_profile.index[_p]] for _p in plot_params])

    * calculating the 1D profile for: param1 with index 0
    * generating samples for profiler
    - number of random search samples = 20000
    - sampling the distribution
    - time taken to sample the distribution 0.09873 (s)
    - doing initial randomized search
    - time taken for random algorithm 0.002526 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 6 / 58
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.244 (s)
      at the end of descent after 100 steps 9 samples were still beeing optimized.
      gradient descents did not improve 39 points
    - doing minimization polishing (scipy)
    - time taken for polishing 7.942 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0005996 (s)
    - smoothing results
    - time taken for smoothing 0.001236 (s)

(3, 1)

Profile accuracy tests:

The profiler calculation is hard. As you might have seen it requires hundreds or thousands of minimization instances.

In this section we investigate the reliability of the profile calculation. We can do so since - in this example - we have the exact distribution available.

We implemented methods to wrap a tensorflow or scipy distribution into a flow so that it can be used for thorough tests.

# In this case we can compare with the profile obtained from the exact distribution:

from tensiometer.synthetic_probability import analytic_flow

exact_flow = analytic_flow.analytic_flow(analytic_flow.tf_prob_wrapper(distribution),
                           param_names=flow.param_names, 
                           param_labels=flow.param_labels, 
                           lims=flow.parameter_ranges)
exact_profile = flow_profiler.posterior_profile_plotter(exact_flow, 
                                          initialize_cache=False,
                                          feedback=2 )
exact_profile.update_cache(params=profile_params, 
                            **profiler_options)

Flow profiler
  * use_scipy  = True
  * pre_polish = True
  * polish     = True
  * smoothing  = True
  * box_prior  = False
  * parameters = ['param1', 'param2', 'param3', 'param4', 'param5', 'param6']
  * feedback   = 2
  * initializing profiler data.
    * generating samples for profiler
    - number of random search samples = 20000
    - sampling the distribution
    - time taken to sample the distribution 0.1024 (s)
  * finding MAP
    * finding global best fit
    - doing initial randomized search
    - time taken for random initial selection 0.0004132 (s)
    - doing minimization (scipy)
    - time taken for minimization 71.38 (s)
  * initializing 1D profiles
    1D profiles:   0%|          | 0/2 [00:00<?, ?it/s]
    * calculating the 1D profile for: param3 with index 2
    - doing initial randomized search
    - time taken for random algorithm 0.0026 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 1 / 63
    - doing gradient descent pre-polishing
    - time taken for gradient descents 18.62 (s)
      at the end of descent after 1 steps 0 samples were still beeing optimized.
      gradient descents did not improve 63 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  47 points
      but  16 samples were still better.
    - time taken for polishing 74.88 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0004458 (s)
    - smoothing results
    - time taken for smoothing 0.0007787 (s)
    1D profiles:  50%|█████     | 1/2 [01:33<01:33, 93.51s/it]
    * calculating the 1D profile for: param5 with index 4
    - doing initial randomized search
    - time taken for random algorithm 0.00265 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 0 / 64
    - doing gradient descent pre-polishing
    - time taken for gradient descents 19 (s)
      at the end of descent after 1 steps 0 samples were still beeing optimized.
      gradient descents did not improve 64 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  44 points
      but  20 samples were still better.
    - time taken for polishing 91.16 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0004303 (s)
    - smoothing results
    - time taken for smoothing 0.0007715 (s)
    1D profiles: 100%|██████████| 2/2 [03:23<00:00, 101.84s/it]
  * initializing 2D profiles
    2D profiles:   0%|          | 0/1 [00:00<?, ?it/s]
    * calculating the 2D profile for: param3, param5 with indexes  2 4
    - doing initial randomized search
    - time taken for random algorithm 0.00413 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 327 / 697
    - doing gradient descent pre-polishing
    - time taken for gradient descents 416.2 (s)
      at the end of descent after 2 steps 0 samples were still beeing optimized.
      gradient descents did not improve 696 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  564 points
      but  133 samples were still better.
    - time taken for polishing 715 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.008266 (s)
    - smoothing results
    - time taken for smoothing 0.00169 (s)
    2D profiles: 100%|██████████| 1/1 [18:51<00:00, 1131.21s/it]

We now plot the profiles.

We log plot the 1D distributions to appreciate flow accuracy in the tails. Note that the flow is optimized on log probabilities so it is expected to do a fairly good job across orders of magnitudes in probability. On the other hand small errors in log space are big errors in probability, hence the requirement of high overall accuracy.

g = plots.get_subplot_plotter()
g.triangle_plot([flow_profile, exact_profile], 
                params=profile_params,
                filled=False, 
                shaded=True, 
                diag1d_kwargs={'normalized':True},
                markers=[exact_profile.flow_MAP[flow_profile.index[_p]] for _p in profile_params],
                legend_labels=['Flow Profile','Exact Profile'])
# log axis on the diagonal:
for _i in range(len(profile_params)):
    _ax = g.subplots[_i, _i]
    _ax.set_yscale('log')
    _ax.set_ylim(1.e-5, 1.e1)
    _ax.set_ylabel('$\\log_{10}(P)$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')

If we trained average flows we can also estimate the error on the profile as the variance of the logPs across the flows. This is what we are going to do here:

g = plots.get_subplot_plotter()
g.triangle_plot([flow_profile, exact_profile], params=profile_params,
                filled=False, shaded=True, diag1d_kwargs={'normalized':True},
                markers=[flow_profile.flow_MAP[flow_profile.index[_p]] for _p in profile_params],
                legend_labels=['Flow Profile','Exact Profile'])

# add error bar on the diagonal:
for _i in range(len(profile_params)):
    # call the method that computes the variance of the profile for an average flow:
    _x, _prob, _temp_std = flow_profile.get_1d_profile_variance(profile_params[_i])
    # do the plotting:
    _ax = g.subplots[_i, _i]
    _ax.plot(_x, _prob, color='k', linestyle='-', label='True')
    _ax.fill_between(_x, _prob - _temp_std, _prob + _temp_std, color='k', alpha=0.2)
    _ax.set_ylabel('$P$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')    

WARNING:tensorflow:5 out of the last 13481 calls to <function FlowCallback.log_probability at 0x7ef9b015a1f0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

WARNING:tensorflow:5 out of the last 13481 calls to <function FlowCallback.log_probability at 0x7ef9b015a1f0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

WARNING:tensorflow:6 out of the last 13482 calls to <function FlowCallback.log_probability at 0x7ef950707ee0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

WARNING:tensorflow:6 out of the last 13482 calls to <function FlowCallback.log_probability at 0x7ef950707ee0> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

Real world application: cosmological parameter profiles

In this section we show a real example of a profile applied to cosmological parameter posteriors.

In this case it is particularly interesting to combine profiling and marginalization. The full parameter space of the model is large - of order 30D - but most of these parameters describe systematic effects. These can be marginalized over, after all we might not really be interested in what happens with them and then we can profile cosmological parameters.

# we start by loading up the posterior:

# load the samples (remove no burn in since the example chains have already been cleaned):
chains_dir = os.path.realpath(os.path.join(os.getcwd(), '../..', 'test_chains'))
# the DES Y1 3x2 chain:
data_chain = getdist.mcsamples.loadMCSamples(file_root=os.path.join(chains_dir, 'DES'), no_cache=True, settings=getdist_settings)

# let's add omegab as a derived parameter:
for _ch in [data_chain]:
    _p = _ch.getParams()
    _h = _p.H0 / 100.
    _ch.addDerived(_p.omegabh2 / _h**2, name='omegab', label='\\Omega_b')
    _ch.updateBaseStatistics()

# we define the parameters of the problem:
param_names = ['H0', 'omegam', 'sigma8', 'ns', 'omegab']

# we then train the flows on the base parameters that we want to use:
kwargs = {
          'feedback': 1,
          'verbose': -1,
          'plot_every': 1000,
          'pop_size': 1,
          'num_flows': 5,
          'epochs': 500,
        }

# actual flow training (note caching):
data_flow = synprob.average_flow_from_chain(data_chain, 
                                            param_names=param_names,
                                            **kwargs)

# plot to make sure training went well:
data_flow.training_plot()

Training flow 0

0epoch [00:00, ?epoch/s]

Training flow 1

0epoch [00:00, ?epoch/s]

Training flow 2

0epoch [00:00, ?epoch/s]

Training flow 3

0epoch [00:00, ?epoch/s]

Training flow 4

0epoch [00:00, ?epoch/s]

# sanity check triangle plot:
g = plots.get_subplot_plotter()
g.triangle_plot([data_chain, data_flow.MCSamples(20000, settings=getdist_settings)], 
                params=param_names,
                filled=False)

# define the options for the profiler:
profiler_options = {
        'num_gd_interactions_1D': 100,  # number of gradient descent interactions for the 1D profile
        'num_gd_interactions_2D': 100,  # number of gradient descent interactions for the 2D profile
        'scipy_options': {  # options for the polishing minimizer
                    'ftol': 1.e-06,
                    'gtol': 0.0,
                    'maxls': 100,
                },
        'scipy_use_jac': True,  # use the jacobian in the minimizer
        'num_points_1D': 64, # number of points for the 1D profile
        'num_points_2D': 32, # number of points per dimension for the 2D profile
        'smooth_scale_1D': 0.2, # smoothing scale for the 1D profile
        'smooth_scale_2D': 0.2, # smoothing scale for the 2D profile
        }

# initialize the profiler:
data_flow_profile = flow_profiler.posterior_profile_plotter(data_flow, 
                                                            initialize_cache=False,  # usually we want to avoid initializing the cache for all the parameters
                                                            feedback=2,  # we want high feedback to see the progress
                                                            )

Flow profiler
  * use_scipy  = True
  * pre_polish = True
  * polish     = True
  * smoothing  = True
  * box_prior  = False
  * parameters = ['H0', 'omegam', 'sigma8', 'ns', 'omegab']
  * feedback   = 2

# now we initialize the cache for the parameters we want to profile:
profile_params = ['omegam', 'sigma8', 'ns']
data_flow_profile.update_cache(params=profile_params, 
                               **profiler_options)

  * initializing profiler data.
    * generating samples for profiler
    - number of random search samples = 20000
    - sampling the distribution
    - time taken to sample the distribution 0.0769 (s)
  * finding MAP
    * finding global best fit
    - doing initial randomized search
    - time taken for random initial selection 0.0005658 (s)
    - doing minimization (scipy)
    - time taken for minimization 9.622 (s)
  * initializing 1D profiles
    1D profiles:   0%|          | 0/3 [00:00<?, ?it/s]
    * calculating the 1D profile for: omegam with index 1
    - doing initial randomized search
    - time taken for random algorithm 0.002404 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 11 / 53
    - doing gradient descent pre-polishing
    - time taken for gradient descents 29.37 (s)
      at the end of descent after 100 steps 42 samples were still beeing optimized.
      gradient descents did not improve 4 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  1 points
      but  52 samples were still better.
    - time taken for polishing 8.643 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0005891 (s)
    - smoothing results
    - time taken for smoothing 0.01072 (s)
    1D profiles:  33%|███▎      | 1/3 [00:38<01:16, 38.04s/it]
    * calculating the 1D profile for: sigma8 with index 2
    - doing initial randomized search
    - time taken for random algorithm 0.002518 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 11 / 53
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.118 (s)
      at the end of descent after 100 steps 45 samples were still beeing optimized.
      gradient descents did not improve 3 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  1 points
      but  52 samples were still better.
    - time taken for polishing 6.939 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.000545 (s)
    - smoothing results
    - time taken for smoothing 0.01069 (s)
    1D profiles:  67%|██████▋   | 2/3 [00:46<00:20, 20.42s/it]
    * calculating the 1D profile for: ns with index 3
    - doing initial randomized search
    - time taken for random algorithm 0.002514 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 0 / 64
    - doing gradient descent pre-polishing
    - time taken for gradient descents 0.8935 (s)
      at the end of descent after 100 steps 59 samples were still beeing optimized.
    - doing minimization polishing (scipy)
    - time taken for polishing 6.413 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0005076 (s)
    - smoothing results
    - time taken for smoothing 0.01098 (s)
    1D profiles: 100%|██████████| 3/3 [00:53<00:00, 17.82s/it]
  * initializing 2D profiles
    2D profiles:   0%|          | 0/3 [00:00<?, ?it/s]
    * calculating the 2D profile for: omegam, sigma8 with indexes  1 2
    - doing initial randomized search
    - time taken for random algorithm 0.004182 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 879 / 145
    - doing gradient descent pre-polishing
    - time taken for gradient descents 0.9864 (s)
      at the end of descent after 100 steps 125 samples were still beeing optimized.
      gradient descents did not improve 6 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  4 points
      but  141 samples were still better.
    - time taken for polishing 19.46 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.002659 (s)
    - smoothing results
    - time taken for smoothing 0.02109 (s)
    2D profiles:  33%|███▎      | 1/3 [00:20<00:41, 20.71s/it]
    * calculating the 2D profile for: omegam, ns with indexes  1 3
    - doing initial randomized search
    - time taken for random algorithm 0.003817 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 465 / 559
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.66 (s)
      at the end of descent after 100 steps 522 samples were still beeing optimized.
      gradient descents did not improve 3 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  1 points
      but  558 samples were still better.
    - time taken for polishing 64.46 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.008546 (s)
    - smoothing results
    - time taken for smoothing 0.02142 (s)
    2D profiles:  67%|██████▋   | 2/3 [01:26<00:47, 47.46s/it]
    * calculating the 2D profile for: sigma8, ns with indexes  2 3
    - doing initial randomized search
    - time taken for random algorithm 0.004856 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 394 / 630
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.896 (s)
      at the end of descent after 100 steps 188 samples were still beeing optimized.
      gradient descents did not improve 2 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  15 points
      but  615 samples were still better.
    - time taken for polishing 39.7 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.009875 (s)
    - smoothing results
    - time taken for smoothing 0.02171 (s)
    2D profiles: 100%|██████████| 3/3 [02:08<00:00, 42.85s/it]

We can now plot the profiles.

Note that - in this case - since we have no evidence estimate available it is crucial to train a bunch of flows to get an estimate of the variance of the profile. If the variance is large it is usually a good idea to retrain...

g = plots.get_subplot_plotter()
g.triangle_plot([data_flow_profile, data_flow.MCSamples(10000)], params=profile_params,
                filled=False, shaded=True, diag1d_kwargs={'normalized':True},
                markers=[data_flow_profile.flow_MAP[data_flow_profile.index[_p]] for _p in profile_params],
                legend_labels=['Profile','Marginal'])

# add error bar on the diagonal:
for _i in range(len(profile_params)):
    _x, _prob, _temp_std = data_flow_profile.get_1d_profile_variance(profile_params[_i])
    # do the plotting:
    _ax = g.subplots[_i, _i]
    _ax.plot(_x, _prob, color='k', linestyle='-', label='True')
    _ax.fill_between(_x, _prob - _temp_std, _prob + _temp_std, color='k', alpha=0.2)
    _ax.set_ylabel('$P$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')    

2D profiles are pretty expensive, 1D profiles, on the other hand, are fairly fast, to the point that we can compute them all for this example:

g = plots.get_subplot_plotter()
data_flow_profile.normalize()
g.plots_1d([data_flow_profile, data_flow.MCSamples(10000)], 
           params=param_names, 
           legend_labels=['Profile','Marginal'],
           markers=[data_flow_profile.flow_MAP[data_flow_profile.index[_p]] for _p in param_names], 
           nx=5, share_y=True, normalize=False)

# add error bars:
for _i in range(len(param_names)):
    _x, _prob, _temp_std = data_flow_profile.get_1d_profile_variance(param_names[_i], normalize_by='max')
    # do the plotting:
    _ax = g.subplots.flatten()[_i]
    _ax.plot(_x, _prob, color='k', linestyle='-', label='True')
    _ax.fill_between(_x, _prob - _temp_std, _prob + _temp_std, color='k', alpha=0.2)
    _ax.set_ylabel('$P$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')

    * calculating the 1D profile for: H0 with index 0
    * generating samples for profiler
    - number of random search samples = 20000
    - sampling the distribution
    - time taken to sample the distribution 0.07819 (s)
    - doing initial randomized search
    - time taken for random algorithm 0.00239 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 8 / 56
    - doing gradient descent pre-polishing
    - time taken for gradient descents 0.9087 (s)
      at the end of descent after 100 steps 56 samples were still beeing optimized.
    - doing minimization polishing (scipy)
    - time taken for polishing 6.728 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0005324 (s)
    - smoothing results
    - time taken for smoothing 0.01068 (s)

    * calculating the 1D profile for: omegab with index 4
    - doing initial randomized search
    - time taken for random algorithm 0.002447 (s)
    - number of 1D bins 64
    - number of empty/filled 1D bins 1 / 63
    - doing gradient descent pre-polishing
    - time taken for gradient descents 0.9839 (s)
      at the end of descent after 100 steps 58 samples were still beeing optimized.
      gradient descents did not improve 4 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  3 points
      but  60 samples were still better.
    - time taken for polishing 10.33 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.0005958 (s)
    - smoothing results
    - time taken for smoothing 0.01128 (s)

If we compute all 1D profiles we can get best-fit and error bars from the profile likelihood (posterior really) ratios. These are defined from posterior thresholds from the maximum.

To do so we have hijacked getdist method getLikeStats.

likestats = data_flow_profile.getLikeStats()
print(likestats)

Best fit sample -log(Like) = -8.509365
Ln(mean 1/like) = -1.220827
mean(-Ln(like)) = 5.682701
-Ln(mean like)  = 10.491468
2*Var(Ln(like)) = 5.342697

parameter   bestfit        lower1         upper1         lower2         upper2
H0          9.8290474E+01  6.7385115E+01  7.5291754E+01  6.1455136E+01  1.0000000E+02   H_0
omegam      2.3610687E-01  1.9742227E-01  2.7118986E-01  1.6668578E-01  3.7569393E-01   \Omega_m
sigma8      9.0459049E-01  7.9663249E-01  9.9866977E-01  6.7035920E-01  1.1249431E+00   \sigma_8
ns          1.0405004E+00  8.3750000E-01  1.2000000E+00  8.0000000E-01  1.2000000E+00   n_s
omegab      7.0236217E-02  5.0739615E-02  8.1745317E-02  1.4963804E-02  9.6055642E-02   \Omega_b

These can be compared to the original margestats. Note that these are obtained starting from the original chain that was fed to the flow:

margestats = data_flow.MCSamples(10000).getMargeStats()
print(margestats)

Marginalized limits: 0.68; 0.95; 0.99

parameter   mean           sddev          lower1         upper1         limit1 lower2         upper2         limit2 lower3         upper3         limit3 
H0          8.4445988E+01  1.2072479E+01  8.0267311E+01  1.0000000E+02  <      6.0797409E+01  1.0000000E+02  <      5.0848907E+01  1.0000000E+02  <       H_0
omegam      2.4949294E-01  3.3890326E-02  2.1366034E-01  2.7202226E-01  two    1.8880543E-01  3.1807320E-01  two    1.8008083E-01  3.7491834E-01  two     \Omega_m
sigma8      8.7898725E-01  7.5195122E-02  8.0734587E-01  9.5208889E-01  two    7.2548592E-01  1.0233500E+00  two    6.6241252E-01  1.0732197E+00  two     \sigma_8
ns          9.9967778E-01  9.8939479E-02  8.8950211E-01  1.1100923E+00  two    8.0000000E-01  1.2000000E+00  none   8.0000000E-01  1.2000000E+00  none    n_s
omegab      5.9108168E-02  1.5951235E-02  4.7156449E-02  7.6493054E-02  two    2.3900382E-02  8.6186484E-02  two    1.5178218E-02  1.0272904E-01  two     \Omega_b

Let's use them to visualize constraints:

g = plots.get_subplot_plotter()
data_flow_profile.normalize()
g.plots_1d([data_flow_profile, data_flow.MCSamples(10000)], 
           params=param_names, legend_labels=['Profile','Marginal'],
           markers=[data_flow_profile.flow_MAP[data_flow_profile.index[_p]] for _p in param_names], 
           nx=5, share_y=True, normalize=False)

# add error bars:
for _i in range(len(param_names)):
    _ax = g.subplots.flatten()[_i]
    _marge = margestats.parWithName(param_names[_i])
    _like = likestats.parWithName(param_names[_i])
    _ax.axvspan(_marge.limits[0].lower, _marge.limits[0].upper, color='r', alpha=0.2)
    _ax.axvspan(_like.ND_limit_bot[0], _like.ND_limit_top[0], color='k', alpha=0.2)

Note that - in this example - confidence intervals obtained from the profile are more conservative than those obained with marginalized statistics.

Having hijacked getdist MCSamples, if we query for MargeStats the profiler we will get confidence intervals calculated with a given mass threshold (i.e. that the isocontour should integrate to a fraction of total).

Let's compare the two:

likestats = data_flow_profile.getLikeStats()
margestats = data_flow_profile.getMargeStats()

g = plots.get_subplot_plotter()
data_flow_profile.normalize()
g.plots_1d([data_flow_profile], 
           params=param_names, legend_labels=['Profile'],
           markers=[data_flow_profile.flow_MAP[data_flow_profile.index[_p]] for _p in param_names], 
           nx=5, share_y=True, normalize=False)

# add error bars:
for _i in range(len(param_names)):
    _ax = g.subplots.flatten()[_i]
    _marge = margestats.parWithName(param_names[_i])
    _like = likestats.parWithName(param_names[_i])
    _ax.axvspan(_marge.limits[0].lower, _marge.limits[0].upper, color='r', alpha=0.2)
    _ax.axvspan(_like.ND_limit_bot[0], _like.ND_limit_top[0], color='k', alpha=0.2)

As we can see the profile threshold is still more conservative. Note that this is a sign of a non-Gaussian distribution and specific to this case.

Best constrained parameters profiles

As discussed in arXiv:2409.09101 projection effects arise because of either true non-Gaussianities of the likelihood or because of unconstrained parameter directions. This means that if we looked at the best constrained directions (that maximize prior to posterior gain, as discussed in arXiv:2112.05737) we have a chance of minimizing projection effects.

If you are interested in how to compute best constrained parameter combinations check out the corresponding notebook.

# let's add S8 as a derived parameter:
for _ch in [data_chain]:
    _p = _ch.getParams()
    _S8 = _p.sigma8 * _p.omegam**0.75
    _ch.addDerived(_S8, name='S8', label='S_8\\equiv \\sigma_8 \\Omega_m^{0.75}')
    _ch.updateBaseStatistics()

# note the slightly different definition of the S8 parameter - that is taken from the notebook

# we define the parameters of the problem:
param_names = ['H0', 'omegam', 'S8', 'ns', 'omegab']

# we then train the flows on the base parameters that we want to use:
kwargs = {
          'feedback': 0,
          'verbose': -1,
          'plot_every': 1000,
          'pop_size': 1,
          'num_flows': 5,
          'epochs': 500,
        }

# actual flow training:
S8_data_flow = synprob.average_flow_from_chain(data_chain, param_names=param_names, **kwargs)

0epoch [00:00, ?epoch/s]

0epoch [00:00, ?epoch/s]

0epoch [00:00, ?epoch/s]

0epoch [00:00, ?epoch/s]

0epoch [00:00, ?epoch/s]

# sanity check triangle plot:
g = plots.get_subplot_plotter()
g.triangle_plot([data_chain, S8_data_flow.MCSamples(20000, settings=getdist_settings),
                 ], 
                params=param_names,
                filled=False)

# define the options for the profiler:
profiler_options = {
        'num_gd_interactions_1D': 100,  # number of gradient descent interactions for the 1D profile
        'num_gd_interactions_2D': 100,  # number of gradient descent interactions for the 2D profile
        'scipy_options': {  # options for the polishing minimizer
                    'ftol': 1.e-06,
                    'gtol': 0.0,
                    'maxls': 100,
                },
        'scipy_use_jac': True,  # use the jacobian in the minimizer
        'num_points_1D': 64, # number of points for the 1D profile
        'num_points_2D': 32, # number of points per dimension for the 2D profile
        'smooth_scale_1D': 0.2, # smoothing scale for the 1D profile
        'smooth_scale_2D': 0.2, # smoothing scale for the 2D profile
        }

# initialize the profiler:
S8_data_flow_profile = flow_profiler.posterior_profile_plotter(S8_data_flow, 
                                                            initialize_cache=False,  # usually we want to avoid initializing the cache for all the parameters
                                                            feedback=1,  # we want high feedback to see the progress
                                                            )

# now we initialize the cache for the parameters we want to profile:
profile_params = ['S8']
S8_data_flow_profile.update_cache(params=profile_params, 
                               **profiler_options)

Flow profiler
  * use_scipy  = True
  * pre_polish = True
  * polish     = True
  * smoothing  = True
  * box_prior  = False
  * parameters = ['H0', 'omegam', 'S8', 'ns', 'omegab']
  * feedback   = 1
  * initializing profiler data.
  * finding MAP
  * initializing 1D profiles
    1D profiles: 100%|██████████| 1/1 [00:37<00:00, 37.78s/it]
  * initializing 2D profiles
    2D profiles: 0it [00:00, ?it/s]

g = plots.get_subplot_plotter()
S8_data_flow_profile.normalize()
g.plots_1d([S8_data_flow_profile, S8_data_flow.MCSamples(10000)], 
           params=profile_params, 
           legend_labels=['Profile','Marginal'],
           markers=[S8_data_flow_profile.flow_MAP[S8_data_flow_profile.index[_p]] for _p in profile_params], 
           nx=5, share_y=True, normalize=False)

# add error bars:
for _i in range(len(profile_params)):
    _x, _prob, _temp_std = S8_data_flow_profile.get_1d_profile_variance(profile_params[_i], normalize_by='max')
    # do the plotting:
    _ax = g.subplots.flatten()[_i]
    _ax.plot(_x, _prob, color='k', linestyle='-', label='True')
    _ax.fill_between(_x, _prob - _temp_std, _prob + _temp_std, color='k', alpha=0.2)
    _ax.set_ylabel('$P$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')

Full profile triangle

By now you probably want to look at a full profile triangle plot...

Let's do it! It might take some time - but appreciate that this was by far impossible to obtain in reasonable times with other methods...

profile_params = data_flow.param_names

g = plots.get_subplot_plotter()
g.triangle_plot([data_flow_profile, data_flow.MCSamples(10000)], params=profile_params,
                filled=False, shaded=True, diag1d_kwargs={'normalized':True},
                markers=[data_flow_profile.flow_MAP[data_flow_profile.index[_p]] for _p in profile_params],
                legend_labels=['Profile','Marginal'])

# add error bar on the diagonal:
for _i in range(len(profile_params)):
    _x, _prob, _temp_std = data_flow_profile.get_1d_profile_variance(profile_params[_i])
    # do the plotting:
    _ax = g.subplots[_i, _i]
    _ax.plot(_x, _prob, color='k', linestyle='-', label='True')
    _ax.fill_between(_x, _prob - _temp_std, _prob + _temp_std, color='k', alpha=0.2)
    _ax.set_ylabel('$P$')
    _ax.tick_params(axis='y', which='both', labelright=True)
    _ax.yaxis.set_label_position('right')    

    * calculating the 2D profile for: H0, omegam with indexes  0 1
    - doing initial randomized search
    - time taken for random algorithm 0.003879 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 579 / 445
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.459 (s)
      at the end of descent after 100 steps 417 samples were still beeing optimized.
      gradient descents did not improve 11 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  5 points
      but  440 samples were still better.
    - time taken for polishing 50.35 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.006629 (s)
    - smoothing results
    - time taken for smoothing 0.02065 (s)

    * calculating the 2D profile for: H0, sigma8 with indexes  0 2
    - doing initial randomized search
    - time taken for random algorithm 0.003905 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 553 / 471
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.514 (s)
      at the end of descent after 100 steps 419 samples were still beeing optimized.
      gradient descents did not improve 9 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  7 points
      but  464 samples were still better.
    - time taken for polishing 53.9 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.006564 (s)
    - smoothing results
    - time taken for smoothing 0.02135 (s)

    * calculating the 2D profile for: H0, ns with indexes  0 3
    - doing initial randomized search
    - time taken for random algorithm 0.004229 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 260 / 764
    - doing gradient descent pre-polishing
    - time taken for gradient descents 2.023 (s)
      at the end of descent after 100 steps 716 samples were still beeing optimized.
    - doing minimization polishing (scipy)
      Warning minimization failed for  4 points
      but  760 samples were still better.
    - time taken for polishing 69.71 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.01042 (s)
    - smoothing results
    - time taken for smoothing 0.02045 (s)

    * calculating the 2D profile for: H0, omegab with indexes  0 4
    - doing initial randomized search
    - time taken for random algorithm 0.004025 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 471 / 553
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.601 (s)
      at the end of descent after 100 steps 545 samples were still beeing optimized.
      gradient descents did not improve 4 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  2 points
      but  551 samples were still better.
    - time taken for polishing 60.38 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.007881 (s)
    - smoothing results
    - time taken for smoothing 0.02098 (s)

    * calculating the 2D profile for: omegam, omegab with indexes  1 4
    - doing initial randomized search
    - time taken for random algorithm 0.003832 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 649 / 375
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.351 (s)
      at the end of descent after 100 steps 52 samples were still beeing optimized.
      gradient descents did not improve 17 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  18 points
      but  357 samples were still better.
    - time taken for polishing 38.12 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.005045 (s)
    - smoothing results
    - time taken for smoothing 0.02113 (s)

    * calculating the 2D profile for: sigma8, omegab with indexes  2 4
    - doing initial randomized search
    - time taken for random algorithm 0.003969 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 605 / 419
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.388 (s)
      at the end of descent after 100 steps 356 samples were still beeing optimized.
      gradient descents did not improve 11 points
    - doing minimization polishing (scipy)
      Warning minimization failed for  12 points
      but  407 samples were still better.
    - time taken for polishing 43.32 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.005517 (s)
    - smoothing results
    - time taken for smoothing 0.02108 (s)

    * calculating the 2D profile for: ns, omegab with indexes  3 4
    - doing initial randomized search
    - time taken for random algorithm 0.004112 (s)
    - number of 2D bins 1024
    - number of empty/filled 2D bins 304 / 720
    - doing gradient descent pre-polishing
    - time taken for gradient descents 1.825 (s)
      at the end of descent after 100 steps 652 samples were still beeing optimized.
    - doing minimization polishing (scipy)
      Warning minimization failed for  3 points
      but  717 samples were still better.
    - time taken for polishing 79.88 (s)
    - doing interpolation on regular grid
    - time taken for interpolation 0.01017 (s)
    - smoothing results
    - time taken for smoothing 0.02086 (s)